#include <stdarg.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "cpgplot.h"

#define BUFSIZE 256

#define RANGE_LOWER 0
#define RANGE_MID   1
#define RANGE_UPPER 2

#define RANGE_TOTAL (RANGE_UPPER+1)

#define RFI_EXTERNAL 0
#define RFI_SELFGEN  1

#define RFI_TOTAL (RFI_SELFGEN+1)

#define COMBINATION_A 0
#define COMBINATION_B 1

#define COMBINATION_TOTAL (COMBINATION_B+1)

#define IF1 0
#define IF2 1

#define IF_TOTAL (IF2+1)

#define POL_A 0
#define POL_B 1
#define POL_COMBINED 2

#define POL_TOTAL (POL_COMBINED+1)

#define ANT1 0
#define ANT2 1
#define ANT3 2
#define ANT4 3
#define ANT5 4
#define ANT6 5

#define ANT_TOTAL (ANT6+1)

#define AVERAGE_SELFGEN_IF1 0
#define AVERAGE_SELFGEN_IF2 1
#define AVERAGE_EXTERNAL_LOW 2
#define AVERAGE_EXTERNAL_HIGH 3

#define AVERAGE_TOTAL (AVERAGE_EXTERNAL_HIGH+1)

#define PLOTOPTIONS_MAKEAXIS 2
#define PLOTOPTIONS_OPENDEVICE 4
#define PLOTOPTIONS_LOGYAXIS 8

#define PLOTTYPE_ONEBASELINE 2
#define PLOTTYPE_ALLBASELINES_SELFGEN 4


typedef float elem_type ;

void set_spectra_colour(int set_or_if_number);
void swap_colour(int *colour_holder);
void plot_setup(char *plot_device,int plot_type,int options);
void topleft_label(char *label,int expected_nlabels,...);
void next_plot(void);
void plot_spectrum(int npoints,float *x_axis,float *y_axis,
		   char *xaxis_title,char *yaxis_title,int options);
elem_type torben(elem_type m[], int n);

int main(int argc,char *argv[]){
  
  float *******frequency;
  float ********amplitude;
  float ******rfi_spectrum;
  float **average_spectra,*cfrequency;
  float tmp_frequency,tmp_amplitude,old_frequency;
  float median_array[3];
  char filename[BUFSIZE],line[BUFSIZE],specref[BUFSIZE],range[BUFSIZE];
  char combination[BUFSIZE],ifn[BUFSIZE],baseline[BUFSIZE],cxaxis[BUFSIZE];
  char pgplot_device[BUFSIZE];
  FILE *read_in=NULL;
  int ******nchans=NULL;
  int nset,wpol,i,j,k,l,m,n,o,plot_option,colour_holder;
  int rfi_mode,curr_combination,curr_if,low_ant,high_ant,cnchans,navged;
  
  /* initialisation */
  nchans=NULL;
  frequency=NULL;
  amplitude=NULL;
  rfi_spectrum=NULL;
  average_spectra=NULL;
  cfrequency=NULL;

  /* put the PGPLOT device you want to use here */
  strcpy(pgplot_device,"cx_rfi.ps/cps");

  /* allocate some memory */
  /* first index is lower antenna */
  nchans=malloc(ANT_TOTAL*sizeof(int*****));
  frequency=malloc(ANT_TOTAL*sizeof(float******));
  amplitude=malloc(ANT_TOTAL*sizeof(float*******));
  rfi_spectrum=malloc(ANT_TOTAL*sizeof(float*****));
  /* second index is higher antenna */
  for (n=ANT1;n<ANT_TOTAL;n++){
    nchans[n]=malloc(ANT_TOTAL*sizeof(int****));
    frequency[n]=malloc(ANT_TOTAL*sizeof(float*****));
    amplitude[n]=malloc(ANT_TOTAL*sizeof(float******));
    rfi_spectrum[n]=malloc(ANT_TOTAL*sizeof(float****));
    /* third index is over all the combinations of IF and frequency */
    for (o=ANT1;o<ANT_TOTAL;o++){
      nchans[n][o]=malloc(COMBINATION_TOTAL*sizeof(int***));
      frequency[n][o]=malloc(COMBINATION_TOTAL*sizeof(float****));
      amplitude[n][o]=malloc(COMBINATION_TOTAL*sizeof(float*****));
      rfi_spectrum[n][o]=malloc(COMBINATION_TOTAL*sizeof(float***));
      /* fourth index is IF */
      for (i=COMBINATION_A;i<COMBINATION_TOTAL;i++){
	nchans[n][o][i]=malloc(IF_TOTAL*sizeof(int**));
	frequency[n][o][i]=malloc(IF_TOTAL*sizeof(float***));
	amplitude[n][o][i]=malloc(IF_TOTAL*sizeof(float****));
	rfi_spectrum[n][o][i]=malloc(IF_TOTAL*sizeof(float**));
	/* fifth index is channel/frequency axis */
	for (j=IF1;j<IF_TOTAL;j++){
	  nchans[n][o][i][j]=malloc(RFI_TOTAL*sizeof(int*));
	  frequency[n][o][i][j]=malloc(RFI_TOTAL*sizeof(float**));
	  amplitude[n][o][i][j]=malloc(RFI_TOTAL*sizeof(float***));
	  rfi_spectrum[n][o][i][j]=malloc(RFI_TOTAL*sizeof(float*));
	  /* sixth index is lower/mid/upper range */
	  for (k=RFI_EXTERNAL;k<RFI_TOTAL;k++){
	    nchans[n][o][i][j][k]=malloc(RANGE_TOTAL*sizeof(int));
	    frequency[n][o][i][j][k]=malloc(RANGE_TOTAL*sizeof(float*));
	    amplitude[n][o][i][j][k]=malloc(RANGE_TOTAL*sizeof(float**));
	    rfi_spectrum[n][o][i][j][k]=NULL;
	    /* various seventh indices */
	    for (l=RANGE_LOWER;l<RANGE_TOTAL;l++){
	      frequency[n][o][i][j][k][l]=NULL; /* channel number */
	      amplitude[n][o][i][j][k][l]=malloc(POL_TOTAL*sizeof(float*)); /* pol */
	      /* eighth index is channel number */
	      for (m=POL_A;m<POL_TOTAL;m++){
		amplitude[n][o][i][j][k][l][m]=NULL;
	      }
	    }
	  }
	}
      }
    }
  }

  /* first we need to loop over the lower antenna */
  for (low_ant=ANT1;low_ant<(ANT_TOTAL-1);low_ant++){
    /* next we loop over the high antenna */
    for (high_ant=low_ant+1;high_ant<ANT_TOTAL;high_ant++){
      /* now we need to loop over all combinations of IF and 
	 frequency */
      for (curr_combination=COMBINATION_A;
	   curr_combination<COMBINATION_TOTAL;curr_combination++){
	switch (curr_combination){
	case COMBINATION_A:
	  strcpy(combination,"a");
	  break;
	case COMBINATION_B:
	  strcpy(combination,"b");
	  break;
	}
	/* now loop over frequency */
	for (curr_if=IF1;curr_if<IF_TOTAL;curr_if++){
	  switch (curr_if){
	  case IF1:
	    strcpy(ifn,"IF1");
	    break;
	  case IF2:
	    strcpy(ifn,"IF2");
	    break;
	  }
	  /* now loop over external and self-generated RFI */
	  for (rfi_mode=0;rfi_mode<RFI_TOTAL;rfi_mode++){
	    switch (rfi_mode){
	    case RFI_EXTERNAL:
	      /* we need to look at the frequency spectra for external
		 RFI */
	      strcpy(specref,"freq");
	      break;
	    case RFI_SELFGEN:
	      /* we need to look at the channel spectra for self-generated
		 RFI */
	      strcpy(specref,"chan");
	      break;
	    }
	    /* now loop over the files that we have */
	    for (nset=0;nset<RANGE_TOTAL;nset++){
	      switch (nset){
	      case RANGE_LOWER:
		strcpy(range,"lower");
		break;
	      case RANGE_MID:
		strcpy(range,"mid");
		break;
	      case RANGE_UPPER:
		strcpy(range,"upper");
		break;
	      }
	      snprintf(filename,BUFSIZE,"%1d%1d_%s_%s.%s.%s.log",low_ant+1,high_ant+1,
		       range,combination,ifn,specref);
	      read_in=fopen(filename,"r");
	      if (read_in==NULL){
		/* couldn't open the file */
		printf("Couldn't open file %s\n",filename);
		/* go to the next possible file */
		continue;
	      }
	      nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]=0;
	      wpol=0;
	      old_frequency=-1.0;
	      while(fgets(line,BUFSIZE,read_in)!=NULL){
		if ((sscanf(line," %f %f",&tmp_frequency,&tmp_amplitude))==2){
		  if (tmp_frequency<old_frequency){
		    wpol++;
		    nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]=0;
		  }
		  nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]++;
		  frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]=
		    realloc(frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset],
			    nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]*
			    sizeof(float));
		  frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]
		    [nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]-1]=tmp_frequency;
		  amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][wpol]=
		    realloc(amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][wpol],
			    nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]*
			    sizeof(float));
		  amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][wpol]
		    [nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]-1]=
		    tmp_amplitude;
		  old_frequency=tmp_frequency;
		}
	      }
	      fclose(read_in);
	  
	      /* next we want to average the two polarisations together to 
		 hopefully improve the signal-to-noise ratio */
	      /* the third polarisation of the amplitude array is 
		 the average of the first two */
	      amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][POL_COMBINED]=
		malloc(nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset]*sizeof(float));
	      for (j=0;j<nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset];j++){
		amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][POL_COMBINED][j]=
		  (amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][POL_A][j]+
		   amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][nset][POL_B][j])/2.0;
	      }
	    }
	    
	    /* OK now we want to get our two sets of RFI. */
	    /* To do this we want to get the median of the three sets. The
	       idea behind this is that if you have external RFI, it will
	       stay at the same frequency, so the median of each frequency
	       channel with external RFI will give a good representation of
	       the power expected from that RFI. However, for RFI that does
	       not stay at the same frequency, but rather the same channel
	       as the frequency shifts, the median of whichever frequencies
	       those features appear should be of the level of the noise,
	       since two out of the three observations should observe noise
	       at that frequency. */
	    /* set up the external RFI array */
	    rfi_spectrum[low_ant][high_ant][curr_combination][curr_if][rfi_mode]=
	      malloc(nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_MID]*
		     sizeof(float));
	    /* we'll only be mapping the frequency range of the mid observations */
	    for (i=0;i<nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_MID];i++){
	      /* we only use the averaged polarisation data */
	      /* make the three element array so we can calculate the median */
	      median_array[0]=
		amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_MID][POL_COMBINED][i];
	      for (j=0;j<nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_LOWER];j++){
		if (frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_LOWER][j]==
		    frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_MID][i]){
		  /* make sure we're getting the same frequency channels */
		  median_array[1]=
		    amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_LOWER][POL_COMBINED][j];
		  break; /* don't need to search anymore */
		}
	      }
	      for (j=0;j<nchans[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_UPPER];j++){
		if (frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_UPPER][j]==
		    frequency[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_MID][i]){
		  /* make sure we're getting the same frequency channels */
		  median_array[2]=
		    amplitude[low_ant][high_ant][curr_combination][curr_if][rfi_mode][RANGE_UPPER][POL_COMBINED][j];
		  break; /* don't need to search anymore */
		}
	      }
	      /* get the median value */
	      rfi_spectrum[low_ant][high_ant][curr_combination][curr_if][rfi_mode][i]=torben(median_array,3);
	    }
	  }
	}
      }
    }
  }
  

  /* make some plots */

  /* plot of external RFI, low frequency */
  plot_setup(pgplot_device,PLOTTYPE_ALLBASELINES_SELFGEN,PLOTOPTIONS_OPENDEVICE);
  for (i=ANT1;i<(ANT_TOTAL-1);i++){
    for (j=i+1;j<ANT_TOTAL;j++){
      next_plot();
      plot_option=PLOTOPTIONS_MAKEAXIS;
      for (k=COMBINATION_A;k<COMBINATION_TOTAL;k++){
	set_spectra_colour(k);
	if (nchans[i][j][k][k][RFI_EXTERNAL][RANGE_MID]>0){
	  plot_spectrum(nchans[i][j][k][k][RFI_EXTERNAL][RANGE_MID],
			frequency[i][j][k][k][RFI_EXTERNAL][RANGE_MID],
			rfi_spectrum[i][j][k][k][RFI_EXTERNAL],
			"Frequency (GHz)","Amplitude (Jy)",plot_option | PLOTOPTIONS_LOGYAXIS);
	  if (plot_option & PLOTOPTIONS_MAKEAXIS){
	    snprintf(baseline,BUFSIZE,"Baseline %d-%d",i+1,j+1);
	    swap_colour(&colour_holder);
	    topleft_label(baseline,0);
	    swap_colour(&colour_holder);
	  }
	  plot_option=0;
	}
      }
    }
  }

  /* plot of external RFI, high frequency */
  plot_setup(pgplot_device,PLOTTYPE_ALLBASELINES_SELFGEN,0);
  for (i=ANT1;i<(ANT_TOTAL-1);i++){
    for (j=i+1;j<ANT_TOTAL;j++){
      next_plot();
      plot_option=PLOTOPTIONS_MAKEAXIS;
      for (k=COMBINATION_A;k<COMBINATION_TOTAL;k++){
	if (k==COMBINATION_A){
	  l=IF2;
	} else if (k==COMBINATION_B){
	  l=IF1;
	}
	set_spectra_colour(l);
	if (nchans[i][j][k][l][RFI_EXTERNAL][RANGE_MID]>0){
	  plot_spectrum(nchans[i][j][k][l][RFI_EXTERNAL][RANGE_MID],
			frequency[i][j][k][l][RFI_EXTERNAL][RANGE_MID],
			rfi_spectrum[i][j][k][l][RFI_EXTERNAL],
			"Frequency (GHz)","Amplitude (Jy)",plot_option | PLOTOPTIONS_LOGYAXIS);
	  if (plot_option & PLOTOPTIONS_MAKEAXIS){
	    snprintf(baseline,BUFSIZE,"Baseline %d-%d",i+1,j+1);
	    swap_colour(&colour_holder);
	    topleft_label(baseline,0);
	    swap_colour(&colour_holder);
	  }
	  plot_option=0;
	}
      }
    }
  }

  /* plot of self-generated RFI, IF1 */
  plot_setup(pgplot_device,PLOTTYPE_ALLBASELINES_SELFGEN,0);
  for (i=ANT1;i<(ANT_TOTAL-1);i++){
    for (j=i+1;j<ANT_TOTAL;j++){
      next_plot();
      plot_option=PLOTOPTIONS_MAKEAXIS;
      for (k=COMBINATION_A;k<COMBINATION_TOTAL;k++){
	set_spectra_colour(k);
	if (nchans[i][j][k][k][RFI_SELFGEN][RANGE_MID]>0){
	  plot_spectrum(nchans[i][j][k][k][RFI_SELFGEN][RANGE_MID],
			frequency[i][j][k][k][RFI_SELFGEN][RANGE_MID],
			rfi_spectrum[i][j][k][k][RFI_SELFGEN],
			"Channel","Amplitude (Jy)",plot_option | PLOTOPTIONS_LOGYAXIS);
	  if (plot_option & PLOTOPTIONS_MAKEAXIS){
	    snprintf(baseline,BUFSIZE,"IF1 Baseline %d-%d",i+1,j+1);
	    swap_colour(&colour_holder);
	    topleft_label(baseline,0);
	    swap_colour(&colour_holder);
	  }
	  plot_option=0;
	}
      }
    }
  }

  /* plot of self-generated RFI, IF2 */
  plot_setup(pgplot_device,PLOTTYPE_ALLBASELINES_SELFGEN,0);
  for (i=ANT1;i<(ANT_TOTAL-1);i++){
    for (j=i+1;j<ANT_TOTAL;j++){
      next_plot();
      plot_option=PLOTOPTIONS_MAKEAXIS;
      for (k=COMBINATION_A;k<COMBINATION_TOTAL;k++){
	if (k==COMBINATION_A){
	  l=IF2;
	} else if (k==COMBINATION_B){
	  l=IF1;
	}
	set_spectra_colour(l);
	if (nchans[i][j][k][l][RFI_SELFGEN][RANGE_MID]>0){
	  plot_spectrum(nchans[i][j][k][l][RFI_SELFGEN][RANGE_MID],
			frequency[i][j][k][l][RFI_SELFGEN][RANGE_MID],
			rfi_spectrum[i][j][k][l][RFI_SELFGEN],
			"Channel","Amplitude (Jy)",plot_option | PLOTOPTIONS_LOGYAXIS);
	  if (plot_option & PLOTOPTIONS_MAKEAXIS){
	    snprintf(baseline,BUFSIZE,"IF2 Baseline %d-%d",i+1,j+1);
	    swap_colour(&colour_holder);
	    topleft_label(baseline,0);
	    swap_colour(&colour_holder);
	  }
	  plot_option=0;
	}
      }
    }
  }

  /* the last plot has four panels, two of which are the self-generated
     RFI in each IF, averaged over antenna and polarisation, and the other
     two are the external RFI in each band, averaged over antenna and
     polarisation */
  plot_setup(pgplot_device,PLOTTYPE_ONEBASELINE,0);
  /* calculate the averages */
  average_spectra=malloc(AVERAGE_TOTAL*sizeof(float*));
  for (i=0;i<AVERAGE_TOTAL;i++){
    /* because the two IFs can have different numbers of channels, we will
       always use IF1 when averaging for the external RFI */
    /* allocate for the number of channels */
    switch (i){
    case AVERAGE_SELFGEN_IF1:
      cnchans=nchans[ANT1][ANT2][COMBINATION_A][IF1][RFI_SELFGEN][RANGE_MID];
      cfrequency=frequency[ANT1][ANT2][COMBINATION_A][IF1][RFI_SELFGEN][RANGE_MID];
      strcpy(cxaxis,"Channel");
      break;
    case AVERAGE_SELFGEN_IF2:
      cnchans=nchans[ANT1][ANT2][COMBINATION_A][IF2][RFI_SELFGEN][RANGE_MID];
      cfrequency=frequency[ANT1][ANT2][COMBINATION_A][IF2][RFI_SELFGEN][RANGE_MID];
      strcpy(cxaxis,"Channel");
      break;
    case AVERAGE_EXTERNAL_LOW:
      cnchans=nchans[ANT1][ANT2][COMBINATION_A][IF1][RFI_EXTERNAL][RANGE_MID];
      cfrequency=frequency[ANT1][ANT2][COMBINATION_A][IF1][RFI_EXTERNAL][RANGE_MID];
      strcpy(cxaxis,"Frequency (GHz)");
      break;
    case AVERAGE_EXTERNAL_HIGH:
      cnchans=nchans[ANT1][ANT2][COMBINATION_B][IF1][RFI_EXTERNAL][RANGE_MID];
      cfrequency=frequency[ANT1][ANT2][COMBINATION_B][IF1][RFI_EXTERNAL][RANGE_MID];
      strcpy(cxaxis,"Frequency (GHz)");
      break;
    }
    average_spectra[i]=calloc(cnchans,sizeof(float));

    /* now do the averaging */
    for (j=0;j<cnchans;j++){
      navged=0;
      for (k=ANT1;k<(ANT_TOTAL-1);k++){
	for (l=k+1;l<ANT_TOTAL;l++){
	  switch (i){
	  case AVERAGE_SELFGEN_IF1:
	    average_spectra[i][j]+=rfi_spectrum[k][l][COMBINATION_A][IF1][RFI_SELFGEN][j];
	    navged++;
	    average_spectra[i][j]+=rfi_spectrum[k][l][COMBINATION_B][IF1][RFI_SELFGEN][j];
	    navged++;
	    break;
	  case AVERAGE_SELFGEN_IF2:
	    average_spectra[i][j]+=rfi_spectrum[k][l][COMBINATION_A][IF2][RFI_SELFGEN][j];
	    navged++;
	    average_spectra[i][j]+=rfi_spectrum[k][l][COMBINATION_B][IF2][RFI_SELFGEN][j];
	    navged++;
	    break;
	  case AVERAGE_EXTERNAL_LOW:
	    average_spectra[i][j]+=rfi_spectrum[k][l][COMBINATION_A][IF1][RFI_EXTERNAL][j];
	    navged++;
	    break;
	  case AVERAGE_EXTERNAL_HIGH:
	    average_spectra[i][j]+=rfi_spectrum[k][l][COMBINATION_B][IF1][RFI_EXTERNAL][j];
	    navged++;
	    break;
	  }
	}
      }
      /* now make the average */
      average_spectra[i][j]/=navged;
    }
    plot_option=PLOTOPTIONS_MAKEAXIS | PLOTOPTIONS_LOGYAXIS;
    next_plot();
    cpgsci(1);
    plot_spectrum(cnchans,cfrequency,average_spectra[i],cxaxis,"Amplitude (Jy)",plot_option);
    switch (i){
    case AVERAGE_SELFGEN_IF1:
      topleft_label("IF1 Self-generated RFI",0);
      break;
    case AVERAGE_SELFGEN_IF2:
      topleft_label("IF2 Self-generated RFI",0);
      break;
    case AVERAGE_EXTERNAL_LOW:
      topleft_label("External RFI",0);
      break;
    case AVERAGE_EXTERNAL_HIGH:
      topleft_label("External RFI",0);
      break;
    }
  }

  cpgclos();
  
  exit(0);
}

void set_spectra_colour(int set_or_if_number){

  switch (set_or_if_number){
  case COMBINATION_A:
    /* same as case IF1: */
    cpgsci(2);
    break;
  case COMBINATION_B:
    /* same as case IF2: */
    cpgsci(4);
    break;
  }
}

void swap_colour(int *colour_holder){
  /* This routine checks if the current colour is 
     the standard foreground colour (1), and if it is
     it sets the colour to the value in 'colour_holder'.
     If the colour is not 1, then this routine puts the
     current colour into 'colour_holder' and sets the
     colour to 1. */
  int currcolour,foreground_colour=1;
  
  cpgqci(&currcolour);
  if (currcolour==foreground_colour){
    /* the colour is the foreground colour */
    cpgsci(*colour_holder);
  } else {
    /* we want to set the colour to the foreground colour */
    *colour_holder=currcolour;
    cpgsci(foreground_colour);
  }

}

void plot_setup(char *plot_device,int plot_type,int options){
  
  if (options & PLOTOPTIONS_OPENDEVICE){
    /* open the PGPLOT device */
    cpgopen(plot_device);
  }

  if (plot_type & PLOTTYPE_ONEBASELINE){
    /* we want 4 panels */
    cpgsubp(2,2);
    /* and big characters */
    cpgsch(1.4);
  } else if (plot_type & PLOTTYPE_ALLBASELINES_SELFGEN){
    /* we want 16 panels */
    cpgsubp(4,4);
    /* and big characters */
    cpgsch(1.8);
  }

  
}
 
void next_plot(void){
  cpgpage();
}

void topleft_label(char *label,int expected_nlabels,...){
  float min_x,max_x,min_y,max_y;
  float string_xlength,string_ylength;
  float curr_xpos,curr_ypos,xinterval,yinterval;
  va_list ap;
  char our_label[BUFSIZE];
  
  /* initialise the argument list */
  va_start(ap,expected_nlabels);

  /* This routine should be called with at least one
     label to be written as 'label'. If more labels
     need to be written, the next argument should be
     the number of labels to be written, excluding the
     'label' label (otherwise pass zero). Then pass
     that many labels to the routine. */

  /* get the axis ranges */
  cpgqwin(&min_x,&max_x,&min_y,&max_y);

  /* get the length of the first label */
  strcpy(our_label,label);
  cpglen(4,our_label,&string_xlength,&string_ylength);

  /* put the label at the top left */
  xinterval=(max_x-min_x)/10;
  yinterval=(max_y-min_y)/10;
  curr_xpos=min_x+xinterval;
  curr_ypos=max_y-yinterval;

  while(1){ /* we'll break out when we're finished */
    /* put the label on the plot */
    cpgptxt(curr_xpos,curr_ypos,0.0,0.0,our_label);
    /* do we have another label? */
    if (expected_nlabels>0){
      strcpy(our_label,va_arg(ap,char*));
      /* move down a line */
      curr_ypos-=yinterval;
      /* get the new string's length */
      cpglen(4,our_label,&string_xlength,&string_ylength);
      /* one less extra label */
      expected_nlabels--;
    } else {
      /* no more labels */
      break;
    }
  }
}

void plot_spectrum(int npoints,float *x_axis,float *y_axis,
		   char *xaxis_title,char *yaxis_title,int options){
  float min_x,max_x,min_y,max_y,*xplot=NULL,*yplot=NULL;
  /* the limits extend parameter multiplies the axis limits
     to give a bit of leeway on the plots */
  float limits_extend_x=1;
  float limits_extend_y=2;
  int i,curr_colour;

  xplot=malloc(npoints*sizeof(float));
  yplot=malloc(npoints*sizeof(float));
  for (i=0;i<npoints;i++){
    xplot[i]=x_axis[i];
    if (options & PLOTOPTIONS_LOGYAXIS)
      yplot[i]=log10f(y_axis[i]);
    else
      yplot[i]=y_axis[i];
  }
  
  if (options & PLOTOPTIONS_MAKEAXIS){
    /* determine the limits */
    min_x=max_x=xplot[0];
    min_y=max_y=yplot[0];
    for (i=1;i<npoints;i++){
      if (xplot[i]<min_x)
	min_x=xplot[i];
      if (xplot[i]>max_x)
	max_x=xplot[i];
      if (yplot[i]<min_y)
	min_y=yplot[i];
      if (yplot[i]>max_y)
	max_y=yplot[i];
    }
    /* extend the limits a bit */
/*     printf("limits before = %f -> %f %f -> %f\n", */
/* 	   min_x,max_x,min_y,max_y); */
    min_x/=limits_extend_x;
    max_x*=limits_extend_x;
/*     if (min_y>0) */
/*       min_y/=limits_extend_y; */
/*     else */
/*       min_y*=limits_extend_y; */
    if (max_y>0)
      max_y*=limits_extend_y;
    else
      max_y/=limits_extend_y;
/*     printf("limits after = %f -> %f %f -> %f\n", */
/* 	   min_x,max_x,min_y,max_y); */

    /* what is the colour index at the moment? */
    cpgqci(&curr_colour);

    /* set up the plot range */
    cpgsci(1);
    cpgswin(min_x,max_x,min_y,max_y);
    if (options & PLOTOPTIONS_LOGYAXIS)
      cpgbox("BCNT",0.0,0,"BCLNT",0.0,0);
    else
      cpgbox("BCNT",0.0,0,"BCNT",0.0,0);

    /* put some axis labels on */
    cpglab(xaxis_title,yaxis_title,"");

    /* return the colour index to what it was when we were
       called */
    cpgsci(curr_colour);
  }

  /* make the plot */
  cpgline(npoints,xplot,yplot);

  /* free our memory */
  free(xplot);
  free(yplot);
}

/*
 * The following code is public domain.
 * Algorithm by Torben Mogensen, implementation by N. Devillard.
 * This code in public domain.
 */

elem_type torben(elem_type m[], int n)
{
    int         i, less, greater, equal;
    elem_type  min, max, guess, maxltguess, mingtguess;

    min = max = m[0] ;
    for (i=1 ; i<n ; i++) {
        if (m[i]<min) min=m[i];
        if (m[i]>max) max=m[i];
    }

    while (1) {
        guess = (min+max)/2;
        less = 0; greater = 0; equal = 0;
        maxltguess = min ;
        mingtguess = max ;
        for (i=0; i<n; i++) {
            if (m[i]<guess) {
                less++;
                if (m[i]>maxltguess) maxltguess = m[i] ;
            } else if (m[i]>guess) {
                greater++;
                if (m[i]<mingtguess) mingtguess = m[i] ;
            } else equal++;
        }
        if (less <= (n+1)/2 && greater <= (n+1)/2) break ; 
        else if (less>greater) max = maxltguess ;
        else min = mingtguess;
    }
    if (less >= (n+1)/2) return maxltguess;
    else if (less+equal >= (n+1)/2) return guess;
    else return mingtguess;
}
